Prediction Modeling
Dani Cosme & Rich Lopez
2019-12-02
load packages
library(tidyverse)source data
source("load_data.R")standardize brain data
# roi distributions
dataset %>%
filter(session == "all") %>%
select(subjectID, con, condition, control, roi, process, meanPE) %>%
unique() %>%
filter(!is.na(roi)) %>%
ggplot(aes(meanPE, fill = roi)) +
geom_density(color = NA) +
facet_wrap(~roi, ncol = 4) +
theme_minimal() +
theme(legend.position = "none")
# dot product distribution
dataset %>%
filter(session == "all") %>%
select(subjectID, con, condition, control, process, test, dotProduct) %>%
unique() %>%
ggplot(aes(dotProduct, fill = test)) +
geom_density(color = NA) +
facet_wrap(~process, ncol = 3, scales = "free") +
theme_minimal()
# standardize
betas_std = betas %>%
group_by(roi, session) %>%
mutate(meanPE_std = scale(meanPE, center = TRUE, scale = TRUE),
meanPE = ifelse(meanPE_std > 3 | meanPE_std < -3, NA, meanPE_std)) %>%
ungroup()
dots_std = dots %>%
group_by(map, test, mask, session) %>%
mutate(dotProduct_std = scale(dotProduct, center = TRUE, scale = TRUE),
dotProduct = ifelse(dotProduct_std > 3 | dotProduct_std < -3, 9999, dotProduct_std)) %>%
ungroup()
# join data frames
dataset = full_join(betas_std, dots_std, by = c("subjectID", "con", "process", "condition", "control", "session")) %>%
mutate(subjectID = as.character(subjectID)) %>%
left_join(., ind_diffs, by = "subjectID") %>%
ungroup() %>%
mutate(subjectID = as.factor(subjectID),
condition = as.factor(condition),
control = as.factor(control),
roi = as.factor(roi),
process = as.factor(process),
test = as.factor(test))tidy data
model_data = dataset %>%
filter(((!grepl("neuralsig", map) & test == "association" & mask == "masked") | (grepl("neuralsig", map) & mask == "unmasked")) &
session == "all" & control == "nature" & condition == "food") %>%
unique() %>%
group_by(process, subjectID) %>%
mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE),
processPE = paste0(process, "PE")) %>%
select(-c(xyz, meanPE, meanPE_std, sdPE, dotProduct, map)) %>%
spread(process, dotProduct_std) %>%
spread(processPE, meanProcessPEstd) %>%
group_by(subjectID) %>%
fill(everything(), .direction = "down") %>%
fill(everything(), .direction = "up") %>%
select(-c(roi, cravingPE, craving_regulationPE, test, age, gender, mask)) %>%
unique()run cross-validated models
BMI
data_ctrl = caret::trainControl(method = "cv", number = 5)
model_bmi_roi = caret::train(bmi ~ rewardPE + valuePE + cognitive_controlPE,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)
model_bmi_corr = caret::train(bmi ~ reward + value + cognitive_control,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)
model_bmi_cr = caret::train(bmi ~ craving_regulation,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)
model_bmi_combined = caret::train(bmi ~ rewardPE + valuePE + cognitive_controlPE + reward + value + cognitive_control + craving_regulation,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)
AIC(model_bmi_roi$finalModel, model_bmi_corr$finalModel, model_bmi_cr$finalModel, model_bmi_combined$finalModel)summary(model_bmi_roi)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.2590 -2.1942 -0.2569 1.6289 14.9842
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 23.0416 0.2557 90.124 <2e-16 ***
## rewardPE -1.0833 0.4651 -2.329 0.0207 *
## valuePE 0.5376 0.3327 1.616 0.1073
## cognitive_controlPE 0.1676 0.4604 0.364 0.7161
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.347 on 248 degrees of freedom
## Multiple R-squared: 0.02536, Adjusted R-squared: 0.01357
## F-statistic: 2.151 on 3 and 248 DF, p-value: 0.09438summary(model_bmi_corr)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.2667 -2.2289 -0.1992 1.5942 15.2710
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.9921 0.2712 84.790 <2e-16 ***
## reward -0.7022 0.8992 -0.781 0.436
## value 0.7088 0.7706 0.920 0.359
## cognitive_control -0.1118 0.4115 -0.272 0.786
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.378 on 248 degrees of freedom
## Multiple R-squared: 0.007218, Adjusted R-squared: -0.004791
## F-statistic: 0.601 on 3 and 248 DF, p-value: 0.6149summary(model_bmi_cr)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.2408 -2.2708 -0.2915 1.6320 15.5800
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.9385 0.2555 89.782 <2e-16 ***
## craving_regulation 0.5522 0.3245 1.702 0.09 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.357 on 250 degrees of freedom
## Multiple R-squared: 0.01145, Adjusted R-squared: 0.007499
## F-statistic: 2.896 on 1 and 250 DF, p-value: 0.09002summary(model_bmi_combined)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.2589 -2.2038 -0.2826 1.8436 14.6980
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.7839 0.3155 72.226 <2e-16 ***
## rewardPE -1.0780 0.5524 -1.951 0.0521 .
## valuePE 0.2161 0.4468 0.484 0.6290
## cognitive_controlPE 1.3243 1.0866 1.219 0.2241
## reward -0.4199 0.8981 -0.468 0.6405
## value 0.7507 0.8080 0.929 0.3538
## cognitive_control -0.8837 0.7803 -1.133 0.2585
## craving_regulation 0.5476 0.3451 1.587 0.1139
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.339 on 244 degrees of freedom
## Multiple R-squared: 0.04526, Adjusted R-squared: 0.01787
## F-statistic: 1.652 on 7 and 244 DF, p-value: 0.1215body fat
data_ctrl = caret::trainControl(method = "cv", number = 5)
model_fat_roi = caret::train(fat ~ rewardPE + valuePE + cognitive_controlPE,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)
model_fat_corr = caret::train(fat ~ reward + value + cognitive_control,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)
model_fat_cr = caret::train(fat ~ craving_regulation,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)
model_fat_combined = caret::train(fat ~ rewardPE + valuePE + cognitive_controlPE + reward + value + cognitive_control + craving_regulation,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)
AIC(model_fat_roi$finalModel, model_fat_corr$finalModel, model_fat_cr$finalModel, model_fat_combined$finalModel)summary(model_fat_roi)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.1941 -5.1750 0.9348 5.3876 19.5811
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.1642 0.6245 40.293 <2e-16 ***
## rewardPE -1.5081 1.1370 -1.326 0.186
## valuePE 0.6589 0.8121 0.811 0.418
## cognitive_controlPE -0.1986 1.1212 -0.177 0.860
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.104 on 245 degrees of freedom
## Multiple R-squared: 0.009479, Adjusted R-squared: -0.00265
## F-statistic: 0.7816 on 3 and 245 DF, p-value: 0.5052summary(model_fat_corr)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.991 -5.246 1.268 5.267 18.942
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.9624 0.6573 37.976 <2e-16 ***
## reward -2.5613 2.1623 -1.185 0.237
## value 1.8902 1.8603 1.016 0.311
## cognitive_control 0.1510 0.9927 0.152 0.879
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.101 on 245 degrees of freedom
## Multiple R-squared: 0.01016, Adjusted R-squared: -0.001957
## F-statistic: 0.8385 on 3 and 245 DF, p-value: 0.4739summary(model_fat_cr)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.532 -5.337 1.246 4.997 19.444
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.3371 0.6214 40.771 <2e-16 ***
## craving_regulation 0.1318 0.7871 0.167 0.867
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.109 on 247 degrees of freedom
## Multiple R-squared: 0.0001135, Adjusted R-squared: -0.003935
## F-statistic: 0.02804 on 1 and 247 DF, p-value: 0.8671summary(model_fat_combined)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.822 -5.097 1.056 5.663 18.591
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.82756 0.77580 32.003 <2e-16 ***
## rewardPE -0.86398 1.35204 -0.639 0.523
## valuePE 0.97665 1.09044 0.896 0.371
## cognitive_controlPE 1.97281 2.65142 0.744 0.458
## reward -2.35581 2.19313 -1.074 0.284
## value 1.52194 1.97758 0.770 0.442
## cognitive_control -0.87829 1.90988 -0.460 0.646
## craving_regulation 0.07419 0.84436 0.088 0.930
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.134 on 241 degrees of freedom
## Multiple R-squared: 0.01823, Adjusted R-squared: -0.01028
## F-statistic: 0.6394 on 7 and 241 DF, p-value: 0.7231SCA
tidy data
dots_sca = dots_std %>%
filter((!grepl("neuralsig", map) & mask == "masked") | (grepl("neuralsig", map) & mask == "unmasked")) %>%
unite(variable, process, test, condition, control, sep = "_", remove = TRUE) %>%
mutate(variable = sprintf("multivariate_%s", variable)) %>%
filter(session == "all") %>%
select(-c(map, con, mask, session, dotProduct)) %>%
left_join(., ind_diffs) %>%
select(-c(sample, DBIC_ID, age, gender))
dots_rest_assoc = dots_sca %>%
filter(grepl("rest", variable) & grepl("snack|meal|dessert|food", variable) & grepl("association", variable)) %>%
spread(variable, dotProduct_std) %>%
unique() %>%
ungroup()
dots_rest_uniform = dots_sca %>%
filter(grepl("rest", variable) & grepl("snack|meal|dessert|food", variable) & grepl("uniformity", variable)) %>%
spread(variable, dotProduct_std) %>%
unique() %>%
ungroup()
dots_nature_assoc = dots_sca %>%
filter(grepl("nature", variable) & grepl("snack|meal|dessert|food", variable) & grepl("association", variable)) %>%
spread(variable, dotProduct_std) %>%
unique() %>%
ungroup()
dots_nature_uniform = dots_sca %>%
filter(grepl("nature", variable) & grepl("snack|meal|dessert|food", variable) & grepl("uniformity", variable)) %>%
spread(variable, dotProduct_std) %>%
unique() %>%
ungroup()
betas_sca = betas_std %>%
group_by(subjectID, process, condition, control) %>%
mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE)) %>%
unite(variable, process, condition, control, sep = "_", remove = TRUE) %>%
mutate(variable = sprintf("univariate_%s", variable)) %>%
filter(session == "all") %>%
select(-c(con, session, xyz, roi, meanPE, sdPE, meanPE_std)) %>%
left_join(., ind_diffs) %>%
select(-c(sample, DBIC_ID, age, gender))
betas_rest = betas_sca %>%
filter(grepl("rest", variable) & grepl("snack|meal|dessert|food", variable)) %>%
unique() %>%
spread(variable, meanProcessPEstd) %>%
ungroup()
betas_nature = betas_sca %>%
filter(grepl("nature", variable) & grepl("snack|meal|dessert|food", variable)) %>%
unique() %>%
spread(variable, meanProcessPEstd) %>%
ungroup()BMI
# set na.action for dredge
options(na.action = "na.fail")betas
if (file.exists("bmi_sca_models_betas.RDS")) {
betas_combined = readRDS("bmi_sca_models_betas.RDS")
} else {
# betas > rest
data = betas_rest %>%
select(-fat) %>%
na.omit()
lm_predictors = paste(names(select(betas_rest, -c(subjectID, bmi, fat))), collapse = " + ")
lm_formula = formula(paste0("bmi ~ ", lm_predictors, collapse = " + "))
full_model = lm(lm_formula,
data = data)
betas_rest_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
# betas > nature
data = betas_nature %>%
select(-fat) %>%
na.omit()
lm_predictors = paste(names(select(betas_nature, -c(subjectID, bmi, fat))), collapse = " + ")
lm_formula = formula(paste0("bmi ~ ", lm_predictors, collapse = " + "))
full_model = lm(lm_formula,
data = data)
betas_nature_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
# combine models
betas_combined = bind_rows(betas_rest_sca, betas_nature_sca) %>%
select(AIC, delta, BIC, df, logLik, weight, `(Intercept)`, everything()) %>%
arrange(AIC)
# save
saveRDS(betas_combined, "bmi_sca_models_betas.RDS")
}# set AIC for null model you want to compare model AIC values to
null = betas_combined %>%
arrange(df) %>%
filter(df == 2)
# tidy for plotting
plot_data = betas_combined %>%
filter(!(df == 2 & delta > .75)) %>%
arrange(AIC) %>%
#slice(1:10000) %>%
mutate(specification = row_number(),
better.fit = ifelse(AIC == null$AIC[1], "equal",
ifelse(AIC < null$AIC[1], "yes", "no"))) %>%
gather(variable, val, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
mutate(multivariate = ifelse(grepl("multivariate", variable) & !is.na(val), 1, NA),
cognitive_control = ifelse(grepl("cognitive_control", variable) & !is.na(val), 1, NA),
reward = ifelse(grepl("reward", variable) & !is.na(val), 1, NA),
value = ifelse(grepl("value", variable) & !is.na(val), 1, NA),
craving = ifelse(grepl("craving_a", variable) & !is.na(val), 1, NA),
craving_regulation = ifelse(grepl("craving_regulation", variable) & !is.na(val), 1, NA),
rest = ifelse(grepl("rest", variable) & !is.na(val), 1, NA),
nature = ifelse(grepl("nature", variable) & !is.na(val), 1, NA),
snack = ifelse(grepl("snack", variable) & !is.na(val), 1, NA),
food = ifelse(grepl("food", variable) & !is.na(val), 1, NA),
meal = ifelse(grepl("meal", variable) & !is.na(val), 1, NA),
dessert = ifelse(grepl("dessert", variable) & !is.na(val), 1, NA),
association_test = ifelse(grepl("association", variable) & !is.na(val), 1, NA),
uniformity_test = ifelse(grepl("uniformity", variable) & !is.na(val), 1, NA)) %>%
select(-c(variable, val)) %>%
unique()
order = plot_data %>%
arrange(AIC) %>%
mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, starts_with("better"))) %>%
filter(!is.na(value)) %>%
group_by(variable) %>%
mutate(order = sum(better.fit.num)) %>%
select(variable, order) %>%
unique()
# variables included in model
variable_names = names(select(plot_data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight))
# plot top panel
a = plot_data %>%
ggplot(aes(specification, AIC, color = better.fit)) +
geom_point(shape = "|", size = 4, alpha = .1) +
geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6",
scale_y_continuous(breaks = scales::pretty_breaks(4)) +
labs(x = "", y = "AIC\n") +
theme_minimal() +
theme(legend.position = "none")
# plot bottom panel
plot_data_b = plot_data %>%
gather(variable, val, eval(variable_names)) %>%
left_join(., order, by = "variable") %>%
mutate(variable = gsub("\\(Intercept\\)", "intercept", variable))
b = plot_data_b %>%
ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
geom_point(data = subset(plot_data_b, !is.na(val)), alpha = .1, size = .25, position = position_jitter(width = 0.25, height = 0.25)) +
scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6",
labs(x = "\nspecification number", y = "variables\n") +
theme_minimal() +
theme(legend.position = "none")
cowplot::plot_grid(a, b, ncol = 1, align = "v")
dots
if (file.exists("bmi_sca_models_dots.RDS")) {
dots_combined = readRDS("bmi_sca_models_dots.RDS")
} else {
# dots association > rest
data = dots_rest_assoc %>%
select(-fat) %>%
na.omit()
lm_predictors = paste(names(select(dots_rest_assoc, -c(subjectID, bmi, fat))), collapse = " + ")
lm_formula = formula(paste0("bmi ~ ", lm_predictors, collapse = " + "))
full_model = lm(lm_formula,
data = data)
clust = parallel::makeCluster(getOption("cl.cores", 7))
invisible(parallel::clusterCall(clust, "library", character.only = TRUE))
parallel::clusterExport(clust, "data")
dots_rest_assoc_sca = MuMIn::pdredge(full_model, cluster = clust, rank = "AIC", extra = "BIC")
parallel::stopCluster(clust)
# dots uniformity > rest
data = dots_rest_uniform %>%
select(-fat) %>%
na.omit()
lm_predictors = paste(names(select(dots_rest_uniform, -c(subjectID, bmi, fat))), collapse = " + ")
lm_formula = formula(paste0("bmi ~ ", lm_predictors, collapse = " + "))
full_model = lm(lm_formula,
data = data)
clust = parallel::makeCluster(getOption("cl.cores", 7))
invisible(parallel::clusterCall(clust, "library", character.only = TRUE))
parallel::clusterExport(clust, "data")
dots_rest_uniform_sca = MuMIn::pdredge(full_model, cluster = clust, rank = "AIC", extra = "BIC")
parallel::stopCluster(clust)
# dots association > nature
data = dots_nature_assoc %>%
select(-fat) %>%
na.omit()
lm_predictors = paste(names(select(dots_nature_assoc, -c(subjectID, bmi, fat))), collapse = " + ")
lm_formula = formula(paste0("bmi ~ ", lm_predictors, collapse = " + "))
full_model = lm(lm_formula,
data = data)
clust = parallel::makeCluster(getOption("cl.cores", 7))
invisible(parallel::clusterCall(clust, "library", character.only = TRUE))
parallel::clusterExport(clust, "data")
dots_nature_assoc_sca = MuMIn::pdredge(full_model, cluster = clust, rank = "AIC", extra = "BIC")
parallel::stopCluster(clust)
# dots uniformity > nature
data = dots_nature_uniform %>%
select(-fat) %>%
na.omit()
lm_predictors = paste(names(select(dots_nature_uniform, -c(subjectID, bmi, fat))), collapse = " + ")
lm_formula = formula(paste0("bmi ~ ", lm_predictors, collapse = " + "))
full_model = lm(lm_formula,
data = data)
clust = parallel::makeCluster(getOption("cl.cores", 7))
invisible(parallel::clusterCall(clust, "library", character.only = TRUE))
parallel::clusterExport(clust, "data")
dots_nature_uniform_sca = MuMIn::pdredge(full_model, cluster = clust, rank = "AIC", extra = "BIC")
parallel::stopCluster(clust)
# combine models
dots_combined = bind_rows(dots_rest_assoc_sca, dots_rest_uniform_sca, dots_nature_assoc_sca, dots_nature_uniform_sca) %>%
select(AIC, delta, BIC, df, logLik, weight, `(Intercept)`, everything()) %>%
arrange(AIC)
# save
saveRDS(dots_combined, "bmi_sca_models_dots.RDS")
}
data_combined = bind_rows(betas_combined, dots_combined) %>%
arrange(AIC)data_combined = readRDS("bmi_sca_models.RDS")combined plot
# set AIC for null model you want to compare model AIC values to
null = data_combined %>%
arrange(df) %>%
filter(df == 2)
# tidy for plotting
plot_data = data_combined %>%
filter(!(df == 2 & delta > .75)) %>%
arrange(AIC) %>%
slice(1:10000) %>%
mutate(specification = row_number(),
better.fit = ifelse(AIC == null$AIC[1], "equal",
ifelse(AIC < null$AIC[1], "yes", "no"))) %>%
gather(variable, val, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
mutate(ROI = ifelse(grepl("univariate", variable) & !is.na(val), 1, NA),
multivariate = ifelse(grepl("multivariate", variable) & !is.na(val), 1, NA),
cognitive_control = ifelse(grepl("cognitive_control", variable) & !is.na(val), 1, NA),
reward = ifelse(grepl("reward", variable) & !is.na(val), 1, NA),
value = ifelse(grepl("value", variable) & !is.na(val), 1, NA),
craving = ifelse(grepl("craving_a", variable) & !is.na(val), 1, NA),
craving_regulation = ifelse(grepl("craving_regulation", variable) & !is.na(val), 1, NA),
rest = ifelse(grepl("rest", variable) & !is.na(val), 1, NA),
nature = ifelse(grepl("nature", variable) & !is.na(val), 1, NA),
snack = ifelse(grepl("snack", variable) & !is.na(val), 1, NA),
food = ifelse(grepl("food", variable) & !is.na(val), 1, NA),
meal = ifelse(grepl("meal", variable) & !is.na(val), 1, NA),
dessert = ifelse(grepl("dessert", variable) & !is.na(val), 1, NA),
association_test = ifelse(grepl("association", variable) & !is.na(val), 1, NA),
uniformity_test = ifelse(grepl("uniformity", variable) & !is.na(val), 1, NA)) %>%
select(-c(variable, val)) %>%
unique()
order = plot_data %>%
arrange(AIC) %>%
mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, starts_with("better"))) %>%
filter(!is.na(value)) %>%
group_by(variable) %>%
mutate(order = sum(better.fit.num)) %>%
select(variable, order) %>%
unique()
# variables included in model
variable_names = names(select(plot_data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight))
# plot top panel
a = plot_data %>%
ggplot(aes(specification, AIC, color = better.fit)) +
geom_point(shape = "|", size = 4, alpha = .1) +
geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
scale_y_continuous(breaks = scales::pretty_breaks(4)) +
labs(x = "", y = "AIC\n") +
theme_minimal() +
theme(legend.position = "none")
# plot bottom panel
b = plot_data %>%
gather(variable, value, eval(variable_names)) %>%
left_join(., order, by = "variable") %>%
mutate(value = ifelse(!is.na(value), "|", ""),
variable = gsub("\\(Intercept\\)", "intercept", variable)) %>%
ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
geom_text(aes(label = value), alpha = .1) +
scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
labs(x = "\nspecification number", y = "variables\n") +
theme_minimal() +
theme(legend.position = "none")
p = cowplot::plot_grid(a, b, ncol = 1, align = "v")
ggsave(p, filename = "~/Desktop/plot.png", width = 25, height = 10)# set AIC for null model you want to compare model AIC values to
null = data_combined %>%
arrange(df) %>%
filter(df == 2)
# tidy for plotting
plot_data = data_combined %>%
filter(!(df == 2 & delta > .75)) %>%
arrange(AIC) %>%
slice(1:10000) %>%
mutate(specification = row_number(),
better.fit = ifelse(AIC == null$AIC[1], "equal",
ifelse(AIC < null$AIC[1], "yes", "no"))) %>%
gather(variable, val, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
mutate(ROI = ifelse(grepl("univariate", variable) & !is.na(val), 1, NA),
multivariate = ifelse(grepl("multivariate", variable) & !is.na(val), 1, NA),
cognitive_control = ifelse(grepl("cognitive_control", variable) & !is.na(val), 1, NA),
reward = ifelse(grepl("reward", variable) & !is.na(val), 1, NA),
value = ifelse(grepl("value", variable) & !is.na(val), 1, NA),
craving = ifelse(grepl("craving_a", variable) & !is.na(val), 1, NA),
craving_regulation = ifelse(grepl("craving_regulation", variable) & !is.na(val), 1, NA),
rest = ifelse(grepl("rest", variable) & !is.na(val), 1, NA),
nature = ifelse(grepl("nature", variable) & !is.na(val), 1, NA),
snack = ifelse(grepl("snack", variable) & !is.na(val), 1, NA),
food = ifelse(grepl("food", variable) & !is.na(val), 1, NA),
meal = ifelse(grepl("meal", variable) & !is.na(val), 1, NA),
dessert = ifelse(grepl("dessert", variable) & !is.na(val), 1, NA),
association_test = ifelse(grepl("association", variable) & !is.na(val), 1, NA),
uniformity_test = ifelse(grepl("uniformity", variable) & !is.na(val), 1, NA)) %>%
select(-c(variable, val)) %>%
unique()
order = plot_data %>%
arrange(AIC) %>%
mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, starts_with("better"))) %>%
filter(!is.na(value)) %>%
group_by(variable) %>%
mutate(order = sum(better.fit.num)) %>%
select(variable, order) %>%
unique()
# variables included in model
variable_names = names(select(plot_data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight))
# plot top panel
a = plot_data %>%
ggplot(aes(specification, AIC, color = better.fit)) +
geom_point(shape = "|", size = 4, alpha = .1) +
geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
scale_y_continuous(breaks = scales::pretty_breaks(4)) +
labs(x = "", y = "AIC\n") +
theme_minimal() +
theme(legend.position = "none")
# plot bottom panel
plot_data_b = plot_data %>%
gather(variable, val, eval(variable_names)) %>%
left_join(., order, by = "variable") %>%
mutate(variable = gsub("\\(Intercept\\)", "intercept", variable))
b = plot_data_b %>%
ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
geom_point(data = subset(plot_data_b, !is.na(val)), alpha = .1, size = .25, position = position_jitter(width = 0.25, height = 0.25)) +
scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6",
labs(x = "\nspecification number", y = "variables\n") +
theme_minimal() +
theme(legend.position = "none")
cowplot::plot_grid(a, b, ncol = 1, align = "v")
summarize number of better fitting models
# plot_data %>%
# gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
# mutate(better.fit = ifelse(better.fit == "yes", 1, 0),
# var.better = ifelse(better.fit == 1 & !is.na(value), 1, 0),
# variable = gsub("\\(Intercept\\)", "intercept", variable),
# variable = gsub("health:", "health x ", variable),
# variable = gsub("health", "food type", variable)) %>%
# group_by(variable) %>%
# summarize(sum.var = sum(var.better, na.rm = TRUE),
# sum.all = sum(better.fit, na.rm = TRUE),
# percent = (sum.var / sum.all) * 100) %>%
# kable(format = "pandoc", digits = 2)% fat
betas
if (file.exists("fat_sca_models_betas.RDS")) {
betas_combined = readRDS("fat_sca_models_betas.RDS")
} else {
# betas > rest
data = betas_rest %>%
select(-bmi) %>%
na.omit()
lm_predictors = paste(names(select(betas_rest, -c(subjectID, bmi, fat))), collapse = " + ")
lm_formula = formula(paste0("fat ~ ", lm_predictors, collapse = " + "))
full_model = lm(lm_formula,
data = data)
betas_rest_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
# betas > nature
data = betas_nature %>%
select(-bmi) %>%
na.omit()
lm_predictors = paste(names(select(betas_nature, -c(subjectID, bmi, fat))), collapse = " + ")
lm_formula = formula(paste0("fat ~ ", lm_predictors, collapse = " + "))
full_model = lm(lm_formula,
data = data)
betas_nature_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
# combine models
betas_combined = bind_rows(betas_rest_sca, betas_nature_sca) %>%
select(AIC, delta, BIC, df, logLik, weight, `(Intercept)`, everything()) %>%
arrange(AIC)
# save
saveRDS(betas_combined, "fat_sca_models_betas.RDS")
}# set AIC for null model you want to compare model AIC values to
null = betas_combined %>%
arrange(df) %>%
filter(df == 2)
# tidy for plotting
plot_data = betas_combined %>%
filter(!(df == 2 & delta > .75)) %>%
arrange(AIC) %>%
#slice(1:10000) %>%
mutate(specification = row_number(),
better.fit = ifelse(AIC == null$AIC[1], "equal",
ifelse(AIC < null$AIC[1], "yes", "no"))) %>%
gather(variable, val, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
mutate(multivariate = ifelse(grepl("multivariate", variable) & !is.na(val), 1, NA),
cognitive_control = ifelse(grepl("cognitive_control", variable) & !is.na(val), 1, NA),
reward = ifelse(grepl("reward", variable) & !is.na(val), 1, NA),
value = ifelse(grepl("value", variable) & !is.na(val), 1, NA),
craving = ifelse(grepl("craving_a", variable) & !is.na(val), 1, NA),
craving_regulation = ifelse(grepl("craving_regulation", variable) & !is.na(val), 1, NA),
rest = ifelse(grepl("rest", variable) & !is.na(val), 1, NA),
nature = ifelse(grepl("nature", variable) & !is.na(val), 1, NA),
snack = ifelse(grepl("snack", variable) & !is.na(val), 1, NA),
food = ifelse(grepl("food", variable) & !is.na(val), 1, NA),
meal = ifelse(grepl("meal", variable) & !is.na(val), 1, NA),
dessert = ifelse(grepl("dessert", variable) & !is.na(val), 1, NA),
association_test = ifelse(grepl("association", variable) & !is.na(val), 1, NA),
uniformity_test = ifelse(grepl("uniformity", variable) & !is.na(val), 1, NA)) %>%
select(-c(variable, val)) %>%
unique()
order = plot_data %>%
arrange(AIC) %>%
mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, starts_with("better"))) %>%
filter(!is.na(value)) %>%
group_by(variable) %>%
mutate(order = sum(better.fit.num)) %>%
select(variable, order) %>%
unique()
# variables included in model
variable_names = names(select(plot_data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight))
# plot top panel
a = plot_data %>%
ggplot(aes(specification, AIC, color = better.fit)) +
geom_point(shape = "|", size = 4, alpha = .1) +
geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6",
scale_y_continuous(breaks = scales::pretty_breaks(4)) +
labs(x = "", y = "AIC\n") +
theme_minimal() +
theme(legend.position = "none")
# plot bottom panel
plot_data_b = plot_data %>%
gather(variable, val, eval(variable_names)) %>%
left_join(., order, by = "variable") %>%
mutate(variable = gsub("\\(Intercept\\)", "intercept", variable))
b = plot_data_b %>%
ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
geom_point(data = subset(plot_data_b, !is.na(val)), alpha = .1, size = .25, position = position_jitter(width = 0.25, height = 0.25)) +
scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6",
labs(x = "\nspecification number", y = "variables\n") +
theme_minimal() +
theme(legend.position = "none")
cowplot::plot_grid(a, b, ncol = 1, align = "v")
dots
if (file.exists("fat_sca_models_dots.RDS")) {
dots_combined = readRDS("fat_sca_models_dots.RDS")
} else {
# dots association > rest
data = dots_rest_assoc %>%
select(-bmi) %>%
na.omit()
lm_predictors = paste(names(select(dots_rest_assoc, -c(subjectID, bmi, fat))), collapse = " + ")
lm_formula = formula(paste0("fat ~ ", lm_predictors, collapse = " + "))
full_model = lm(lm_formula,
data = data)
dots_rest_assoc_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
# dots uniformity > rest
data = dots_rest_uniform %>%
select(-bmi) %>%
na.omit()
lm_predictors = paste(names(select(dots_rest_uniform, -c(subjectID, bmi, fat))), collapse = " + ")
lm_formula = formula(paste0("fat ~ ", lm_predictors, collapse = " + "))
full_model = lm(lm_formula,
data = data)
dots_rest_uniform_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
# dots association > nature
data = dots_nature_assoc %>%
select(-bmi) %>%
na.omit()
lm_predictors = paste(names(select(dots_nature_assoc, -c(subjectID, bmi, fat))), collapse = " + ")
lm_formula = formula(paste0("fat ~ ", lm_predictors, collapse = " + "))
full_model = lm(lm_formula,
data = data)
dots_nature_assoc_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
# dots association > nature
data = dots_nature_uniform %>%
select(-bmi) %>%
na.omit()
lm_predictors = paste(names(select(dots_nature_uniform, -c(subjectID, bmi, fat))), collapse = " + ")
lm_formula = formula(paste0("fat ~ ", lm_predictors, collapse = " + "))
full_model = lm(lm_formula,
data = data)
dots_nature_uniform_sca = MuMIn::dredge(full_model, rank = "AIC", extra = "BIC")
# combine models
dots_combined = bind_rows(dots_rest_assoc_sca, dots_rest_uniform_sca, dots_nature_assoc_sca, dots_nature_uniform_sca) %>%
select(AIC, delta, BIC, df, logLik, weight, `(Intercept)`, everything()) %>%
arrange(AIC)
# save
saveRDS(dots_combined, "fat_sca_models_dots.RDS")
}
data_combined = bind_rows(betas_combined, dots_combined) %>%
arrange(AIC)# set AIC for null model you want to compare model AIC values to
null = dots_combined %>%
arrange(df) %>%
filter(df == 2)
# tidy for plotting
plot_data = dots_combined %>%
filter(!(df == 2 & delta > .75)) %>%
arrange(AIC) %>%
slice(1:10000) %>%
mutate(specification = row_number(),
better.fit = ifelse(AIC == null$AIC[1], "equal",
ifelse(AIC < null$AIC[1], "yes", "no"))) %>%
gather(variable, val, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
mutate(ROI = ifelse(grepl("univariate", variable) & !is.na(val), 1, NA),
multivariate = ifelse(grepl("multivariate", variable) & !is.na(val), 1, NA),
cognitive_control = ifelse(grepl("cognitive_control", variable) & !is.na(val), 1, NA),
reward = ifelse(grepl("reward", variable) & !is.na(val), 1, NA),
value = ifelse(grepl("value", variable) & !is.na(val), 1, NA),
craving = ifelse(grepl("craving_a", variable) & !is.na(val), 1, NA),
craving_regulation = ifelse(grepl("craving_regulation", variable) & !is.na(val), 1, NA),
rest = ifelse(grepl("rest", variable) & !is.na(val), 1, NA),
nature = ifelse(grepl("nature", variable) & !is.na(val), 1, NA),
snack = ifelse(grepl("snack", variable) & !is.na(val), 1, NA),
food = ifelse(grepl("food", variable) & !is.na(val), 1, NA),
meal = ifelse(grepl("meal", variable) & !is.na(val), 1, NA),
dessert = ifelse(grepl("dessert", variable) & !is.na(val), 1, NA),
association_test = ifelse(grepl("association", variable) & !is.na(val), 1, NA),
uniformity_test = ifelse(grepl("uniformity", variable) & !is.na(val), 1, NA)) %>%
select(-c(variable, val)) %>%
unique()
order = plot_data %>%
arrange(AIC) %>%
mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, starts_with("better"))) %>%
filter(!is.na(value)) %>%
group_by(variable) %>%
mutate(order = sum(better.fit.num)) %>%
select(variable, order) %>%
unique()
# variables included in model
variable_names = names(select(plot_data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight))
# plot top panel
a = plot_data %>%
ggplot(aes(specification, AIC, color = better.fit)) +
geom_point(shape = "|", size = 4, alpha = .1) +
geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
scale_y_continuous(breaks = scales::pretty_breaks(4)) +
labs(x = "", y = "AIC\n") +
theme_minimal() +
theme(legend.position = "none")
# plot bottom panel
plot_data_b = plot_data %>%
gather(variable, val, eval(variable_names)) %>%
left_join(., order, by = "variable") %>%
mutate(variable = gsub("\\(Intercept\\)", "intercept", variable))
b = plot_data_b %>%
ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
geom_point(data = subset(plot_data_b, !is.na(val)), alpha = .1, size = .25, position = position_jitter(width = 0.25, height = 0.25)) +
scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6",
labs(x = "\nspecification number", y = "variables\n") +
theme_minimal() +
theme(legend.position = "none")
cowplot::plot_grid(a, b, ncol = 1, align = "v")combined plot
# set AIC for null model you want to compare model AIC values to
null = data_combined %>%
arrange(df) %>%
filter(df == 2)
# tidy for plotting
plot_data = data_combined %>%
filter(!(df == 2 & delta > .75)) %>%
arrange(AIC) %>%
slice(1:10000) %>%
mutate(specification = row_number(),
better.fit = ifelse(AIC == null$AIC[1], "equal",
ifelse(AIC < null$AIC[1], "yes", "no"))) %>%
gather(variable, val, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
mutate(ROI = ifelse(grepl("univariate", variable) & !is.na(val), 1, NA),
multivariate = ifelse(grepl("multivariate", variable) & !is.na(val), 1, NA),
cognitive_control = ifelse(grepl("cognitive_control", variable) & !is.na(val), 1, NA),
reward = ifelse(grepl("reward", variable) & !is.na(val), 1, NA),
value = ifelse(grepl("value", variable) & !is.na(val), 1, NA),
craving = ifelse(grepl("craving_a", variable) & !is.na(val), 1, NA),
craving_regulation = ifelse(grepl("craving_regulation", variable) & !is.na(val), 1, NA),
rest = ifelse(grepl("rest", variable) & !is.na(val), 1, NA),
nature = ifelse(grepl("nature", variable) & !is.na(val), 1, NA),
snack = ifelse(grepl("snack", variable) & !is.na(val), 1, NA),
food = ifelse(grepl("food", variable) & !is.na(val), 1, NA),
meal = ifelse(grepl("meal", variable) & !is.na(val), 1, NA),
dessert = ifelse(grepl("dessert", variable) & !is.na(val), 1, NA),
association_test = ifelse(grepl("association", variable) & !is.na(val), 1, NA),
uniformity_test = ifelse(grepl("uniformity", variable) & !is.na(val), 1, NA)) %>%
select(-c(variable, val)) %>%
unique()
order = plot_data %>%
arrange(AIC) %>%
mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, starts_with("better"))) %>%
filter(!is.na(value)) %>%
group_by(variable) %>%
mutate(order = sum(better.fit.num)) %>%
select(variable, order) %>%
unique()
# variables included in model
variable_names = names(select(plot_data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight))
# plot top panel
a = plot_data %>%
ggplot(aes(specification, AIC, color = better.fit)) +
geom_point(shape = "|", size = 4, alpha = .1) +
geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
scale_y_continuous(breaks = scales::pretty_breaks(4)) +
labs(x = "", y = "AIC\n") +
theme_minimal() +
theme(legend.position = "none")
# plot bottom panel
plot_data_b = plot_data %>%
gather(variable, val, eval(variable_names)) %>%
left_join(., order, by = "variable") %>%
mutate(variable = gsub("\\(Intercept\\)", "intercept", variable))
b = plot_data_b %>%
ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
geom_point(data = subset(plot_data_b, !is.na(val)), alpha = .1, size = .25, position = position_jitter(width = 0.25, height = 0.25)) +
scale_color_manual(values = c("black", "#F43C13")) + #"#5BBCD6",
labs(x = "\nspecification number", y = "variables\n") +
theme_minimal() +
theme(legend.position = "none")
cowplot::plot_grid(a, b, ncol = 1, align = "v")summarize number of better fitting models
# plot_data %>%
# gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, specification, better.fit)) %>%
# mutate(better.fit = ifelse(better.fit == "yes", 1, 0),
# var.better = ifelse(better.fit == 1 & !is.na(value), 1, 0),
# variable = gsub("\\(Intercept\\)", "intercept", variable),
# variable = gsub("health:", "health x ", variable),
# variable = gsub("health", "food type", variable)) %>%
# group_by(variable) %>%
# summarize(sum.var = sum(var.better, na.rm = TRUE),
# sum.all = sum(better.fit, na.rm = TRUE),
# percent = (sum.var / sum.all) * 100) %>%
# kable(format = "pandoc", digits = 2)